traits
above- and belowground
acquisitive vs. conservative resource-use strategies
collaborative vs. individual resource-use strategies
fine root proportion
thickness of roots
Goal: Examine how traits shift in community compositional (ordinational) space, placing each trait on an axis.
Author: Caryn D. Iwanaga
Updated: 02/06/2025
library(tidyverse)
library(dplyr)
library(ggplot2)
library(httr) # read out Dropbox folders
library(vegan)
library(readxl)
cover.dat.labels <- read.csv("https://www.dropbox.com/scl/fi/up8nnzkcpchsr45f8cm92/Compost_Cover_LongClean.csv?rlkey=z2tvaj8t6khadef7ydz782zka&st=qwef9ys0&dl=1") %>%
mutate(
nut_trt = factor(ifelse(nut_trt =="C", "+compost", ifelse(nut_trt =="F", "+N fertilizer", ifelse(nut_trt =="N", "control", NA))), levels = c("control", "+N fertilizer", "+compost")),
ppt_trt = ifelse(ppt_trt == "D", "dry", ifelse(ppt_trt == "XC", "control", ifelse(ppt_trt == "W", "wet", NA))))
# Define a reusable color scale
# custom_colors <- scale_fill_manual(
# values = c(
# "control" = rgb(195, 197, 193, maxColorValue = 250),
# "+N fertilizer" = rgb(87, 62, 92, maxColorValue = 225),
# "+compost" = rgb(167, 156, 109, maxColorValue = 225)
# )
# )
custom_colors <- scale_fill_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
# Set a global theme
theme_set(theme_bw())
site.df <-
cover.dat.labels[, c(1:7, 18:19)] %>% #subset columns
mutate(
grazing_hist = ifelse(block == 1 | block == 2, "high", ifelse(block == 3 | block == 4, "low", NA))
) %>% # separate by block (grazing intensities)
group_by(nut_trt, ppt_trt, yr, code4, grazing_hist) %>% # make each code unique -- one value for species abundance for each trt
summarise(
# pct_cover,
abundance = mean(pct_cover)
# sd = sd(pct_cover)
)
`summarise()` has grouped output by 'nut_trt', 'ppt_trt', 'yr', 'code4'. You can override using the `.groups` argument.
matrix0 <- spread(
site.df,
code4,
abundance,
fill = 0)
# make rownames ----
rownames(matrix0) <- paste(matrix0$nut_trt, matrix0$ppt_trt, matrix0$yr, matrix0$grazing_hist, sep = "_")
Warning: Setting row names on a tibble is deprecated.
# delete columns ----
matrix <- matrix0[, c(5:79)]
# relativize data with Bray-Curtis dissimilarity
matrix.bray <- decostand(matrix, "total")
rownames(matrix.bray) <- rownames(matrix0)
site.block.df <-
cover.dat.labels[, c(1:7, 18:19)] %>% #subset columns
group_by(nut_trt, ppt_trt, yr, code4, block) %>% # make each code unique -- one value for species abundance for each trt
summarise(
# pct_cover,
abundance = mean(pct_cover)
# sd = sd(pct_cover)
)
`summarise()` has grouped output by 'nut_trt', 'ppt_trt', 'yr', 'code4'. You can override using the `.groups` argument.
site.block.sd.df <-
cover.dat.labels[, c(1:7, 18:19)] %>% #subset columns
group_by(nut_trt, ppt_trt, yr, code4, block) %>% # make each code unique -- one value for species abundance for each trt
summarise(
pct_cover,
abundance = mean(pct_cover),
sd = sd(pct_cover)
)
Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in dplyr 1.1.0.
Please use `reframe()` instead.
When switching from `summarise()` to `reframe()`, remember that `reframe()` always returns an ungrouped data frame and adjust accordingly.`summarise()` has grouped output by 'nut_trt', 'ppt_trt', 'yr', 'code4', 'block'. You can override using the `.groups` argument.
matrix0.block <- spread(
site.block.df,
code4,
abundance,
fill = 0)
# make rownames ----
rownames(matrix0.block) <- paste(matrix0.block$nut_trt, matrix0.block$ppt_trt, matrix0.block$yr, matrix0.block$block, sep = "_")
Warning: Setting row names on a tibble is deprecated.
# delete columns ----
matrix.block <- matrix0.block[, -c(0:4)]
# relativize data with Bray-Curtis dissimilarity
matrix.block.bray <- decostand(matrix.block, "total")
rownames(matrix.block.bray) <- rownames(matrix0.block)
Rules of thumb:
< 0.05 is excellent
0.05-0.1 is good
0.1-0.2 is fair
0.2-0.3 is cause for concern…
> 0.3 is poor, random
Stress: fair
0.1966553
Grouping by year and grazing history had most significant clustering shown.
nmds.comp <- metaMDS(matrix.bray, k=2, trymax = 25)
Run 0 stress 0.2020736
Run 1 stress 0.2044105
Run 2 stress 0.2084765
Run 3 stress 0.2084764
Run 4 stress 0.2026778
Run 5 stress 0.2018338
... New best solution
... Procrustes: rmse 0.06193884 max resid 0.2323138
Run 6 stress 0.1971149
... New best solution
... Procrustes: rmse 0.02534655 max resid 0.1298511
Run 7 stress 0.2027879
Run 8 stress 0.20163
Run 9 stress 0.2033608
Run 10 stress 0.2017485
Run 11 stress 0.2176904
Run 12 stress 0.2013629
Run 13 stress 0.1966668
... New best solution
... Procrustes: rmse 0.01514645 max resid 0.1046137
Run 14 stress 0.1991103
Run 15 stress 0.202025
Run 16 stress 0.2033208
Run 17 stress 0.2020741
Run 18 stress 0.2034839
Run 19 stress 0.2043941
Run 20 stress 0.1991947
Run 21 stress 0.2128964
Run 22 stress 0.1991946
Run 23 stress 0.2039391
Run 24 stress 0.2133823
Run 25 stress 0.196656
... New best solution
... Procrustes: rmse 0.004511995 max resid 0.02760931
*** Best solution was not repeated -- monoMDS stopping criteria:
24: stress ratio > sratmax
1: scale factor of the gradient < sfgrmin
nmds.comp # pretty not great stress
Call:
metaMDS(comm = matrix.bray, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray
Distance: bray
Dimensions: 2
Stress: 0.196656
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 25 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray’
stressplot(nmds.comp)
plot(nmds.comp, type="t")
nmds1 <- as.data.frame(scores(nmds.comp, choices=c(1), display=c("sites"))) %>%
mutate(matrix0[, c(0:4)])
nmds2 <- as.data.frame(scores(nmds.comp, choices=c(2), display=c("sites")))
nmds_dat <- cbind(nmds1, nmds2) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_dat, aes(NMDS1, NMDS2, fill = nut_trt, color = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Nutrient Treatment") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
Please use `linewidth` instead.
ggplot(nmds_dat, aes(NMDS1, NMDS2, color = as.factor(yr), fill = as.factor(yr))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Year")
ggplot(nmds_dat, aes(NMDS1, NMDS2, fill = as.factor(ppt_trt), color = as.factor(ppt_trt))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Precipitation Treatment") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
ggplot(nmds_dat, aes(NMDS1, NMDS2, fill = as.factor(grazing_hist), color = as.factor(grazing_hist))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Grazing History")
ggplot(nmds_dat, aes(NMDS1, NMDS2, color = as.factor(grazing_hist), fill = as.factor(yr))) +
geom_point(aes(fill = as.factor(yr)), shape = 21, size = 4, stroke = 1) +
stat_ellipse(aes(color = as.factor(yr)),
geom = "polygon",
fill = NA,
size = 1) +
labs(title = "Grouped by Grazing History & Year") +
scale_color_manual(
values = c(
"low" = "cyan",
"high" = "brown",
"2019" = "salmon",
"2020" = "limegreen",
"2021" = "cornflowerblue")
)
ggplot(nmds_dat, aes(NMDS1, NMDS2, color = nut_trt, fill = as.factor(yr))) +
geom_point(shape = 21, size = 4, stroke = 1) +
stat_ellipse(aes(color = as.factor(nut_trt)),
geom = "polygon",
fill = NA, # Makes the ellipse transparent inside
size = 1) +
labs(title = "Grouped by Nutrient Treatment & Year") +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
Stress: cause for concern
0.2329306
nmds.block.comp <- metaMDS(matrix.block.bray, k=2, trymax = 25)
nmds.block.comp
stressplot(nmds.block.comp)
plot(nmds.block.comp, type="t")
nmds1 <- as.data.frame(scores(nmds.block.comp, choices=c(1), display=c("sites"))) %>%
mutate(matrix0.block[, c(0:4)])
nmds2 <- as.data.frame(scores(nmds.block.comp, choices=c(2), display=c("sites")))
nmds_dat.block <- cbind(nmds1, nmds2) %>%
select(nut_trt:block, NMDS1, NMDS2)
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, fill = nut_trt, color = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Nutrient Treatment - Blocks") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, color = as.factor(yr), fill = as.factor(yr))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Year - Blocks")
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, fill = as.factor(ppt_trt), color = as.factor(ppt_trt))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Precipitation Treatment - Blocks") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, fill = as.factor(block), color = as.factor(block))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Block")
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, color = as.factor(block), fill = as.factor(yr))) +
geom_point(aes(fill = as.factor(yr)), shape = 21, size = 4, stroke = 1) +
stat_ellipse(aes(color = as.factor(block)),
geom = "polygon",
fill = NA,
size = 1) +
labs(title = "Grouped by Block & Year - Blocks") +
scale_color_manual(
values = c(
# "low" = "cyan",
# "high" = "brown",
"2019" = "salmon",
"2020" = "limegreen",
"2021" = "cornflowerblue")
)
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, color = nut_trt, fill = as.factor(yr))) +
geom_point(shape = 21, size = 4, stroke = 1) +
stat_ellipse(aes(color = as.factor(nut_trt)),
geom = "polygon",
fill = NA, # Makes the ellipse transparent inside
size = 1) +
labs(title = "Grouped by Nutrient Treatment & Year - Blocks") +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
Stress: fair
0.1869068
# relativize data with Bray-Curtis dissimilarity
matrix.low <- matrix0 %>%
filter(grazing_hist == "low")
row.low <- paste(matrix.low$nut_trt, matrix.low$ppt_trt, matrix.low$yr, matrix.low$grazing_hist, sep = "_")
matrix.low <- matrix.low[, c(5:79)]
rownames(matrix.low) <- row.low
matrix.bray.low <- decostand(matrix.low, "total")
nmds.comp.low <- metaMDS(matrix.bray.low, k=2, trymax = 25)
nmds.comp.low # pretty not great stress
stressplot(nmds.comp.low)
plot(nmds.comp.low, type="t")
nmds1 <- as.data.frame(scores(nmds.comp.low, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.low), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.low), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.low), "_"), `[`, 3),
grazing_hist = sapply(strsplit(rownames(matrix.low), "_"), `[`, 4)
)
nmds2 <- as.data.frame(scores(nmds.comp.low, choices=c(2), display=c("sites")))
nmds_low_dat <- cbind(nmds1, nmds2) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_low_dat, aes(NMDS1, NMDS2, color = as.factor(yr), shape = nut_trt)) +
geom_point(size = 3, stroke = 1) +
labs(title = "2019-2021 (LOW) - Grouping by Year & Nutrient Treatment")
ggplot(nmds_low_dat, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
labs(title = "2019-2021 (LOW) - Grouping by Nutrient Treatment")+
custom_colors
ggplot(nmds_low_dat, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2019-2021 (LOW) - Grouping by Precipitation Treatment")
# add year ellipses ----
ggplot(nmds_low_dat, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(aes(color = as.factor(yr)),
geom = "polygon",
fill = NA,
size = 1) +
labs(title = "2019-2021 (LOW) - Grouping by Nutrient Treatment")+
custom_colors
ggplot(nmds_low_dat, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(aes(color = as.factor(yr)),
geom = "polygon",
fill = NA,
size = 1) +
custom_colors +
labs(title = "2019-2021 (LOW) - Grouping by Precipitation Treatment")
# relativize data with Bray-Curtis dissimilarity
matrix.2019 <- matrix0 %>%
filter(yr == 2019)
row.2019 <- paste(matrix.2019$nut_trt, matrix.2019$ppt_trt, matrix.2019$yr, matrix.2019$grazing_hist, sep = "_")
matrix.2019 <- matrix.2019[, c(5:79)]
rownames(matrix.2019) <- row.2019
matrix.bray.2019 <- decostand(matrix.2019, "total")
Stress: Excellent
0.0560178
nmds.2019 <- metaMDS(matrix.bray.2019, k=2, trymax = 25)
nmds.2019 # pretty not great stress
stressplot(nmds.2019)
plot(nmds.2019, type="t")
nmds1.2019 <- as.data.frame(scores(nmds.2019, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 3),
grazing_hist = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 4)
)
nmds2.2019 <- as.data.frame(scores(nmds.2019, choices=c(2), display=c("sites")))
nmds_2019_dat <- cbind(nmds1.2019, nmds2.2019) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_2019_dat, aes(NMDS1, NMDS2, color = grazing_hist, shape = nut_trt)) +
geom_point(size = 3, stroke = 1) +
labs(title = "2019 - Grouping by Grazing History & Nutrient Treatment")
ggplot(nmds_2019_dat, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
labs(title = "2019 - Grouping by Nutrient Treatment")+
custom_colors
ggplot(nmds_2019_dat, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2019 - Grouping by Precipitation Treatment")
# relativize data with Bray-Curtis dissimilarity
matrix.2020 <- matrix0 %>%
filter(yr == 2020)
row.2020 <- paste(matrix.2020$nut_trt, matrix.2020$ppt_trt, matrix.2020$yr, matrix.2020$grazing_hist, sep = "_")
matrix.2020 <- matrix.2020[, c(5:79)]
rownames(matrix.2020) <- row.2020
matrix.bray.2020 <- decostand(matrix.2020, "total")
Stress: fair
0.1522495
nmds.2020 <- metaMDS(matrix.bray.2020, k=2, trymax = 25)
nmds.2020
stressplot(nmds.2020)
plot(nmds.2020, type="t")
nmds1.2020 <- as.data.frame(scores(nmds.2020, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 3),
grazing_hist = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 4)
)
nmds2.2020 <- as.data.frame(scores(nmds.2020, choices=c(2), display=c("sites")))
nmds_2020_dat <- cbind(nmds1.2020, nmds2.2020) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_2020_dat, aes(NMDS1, NMDS2, color = grazing_hist, shape = nut_trt)) +
geom_point(size = 4, stroke = 1) +
labs(title = "2020 - Grouping by Grazing History & Nutrient Treatment")
ggplot(nmds_2020_dat, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2020 - Grouping by Nutrient Treatment")
ggplot(nmds_2020_dat, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2020 - Grouping by Precipitation Treatment")
# relativize data with Bray-Curtis dissimilarity
matrix.2021 <- matrix0 %>%
filter(yr == 2021)
row.2021 <- paste(matrix.2021$nut_trt, matrix.2021$ppt_trt, matrix.2021$yr, matrix.2021$grazing_hist, sep = "_")
matrix.2021 <- matrix.2021[, c(5:79)]
rownames(matrix.2021) <- row.2021
matrix.bray.2021 <- decostand(matrix.2021, "total")
Stress
0.08067789
# relativize data with Bray-Curtis dissimilarity
matrix.2021.1 <- matrix0 %>%
filter(yr == 2021,
grazing_hist == "low")
row.2021.1 <- paste(matrix.2021.1$nut_trt, matrix.2021.1$ppt_trt, matrix.2021.1$yr, matrix.2021.1$grazing_hist, sep = "_")
matrix.2021.1 <- matrix.2021.1[, c(5:79)]
rownames(matrix.2021.1) <- row.2021.1
matrix.bray.2021.1 <- decostand(matrix.2021.1, "total")
nmds.2021.1 <- metaMDS(matrix.bray.2021.1, k=2, trymax = 25)
nmds.2021.1
stressplot(nmds.2021.1)
plot(nmds.2021.1, type="t")
dummy <- as.data.frame(scores(nmds.2021.1, choices=c(1), display=c("sites")))
nmds1.2021.1.df <- as.data.frame(scores(nmds.2021.1, choices=c(1), display=c("sites"))) %>%
mutate(trt = rownames(dummy)) %>%
separate(trt, c("nut_trt", "ppt_trt", "pr", "grazing_hist"), sep="_")
nmds2.2021.1 <- as.data.frame(scores(nmds.2021.1, choices=c(2), display=c("sites")))
nmds_2021_dat.1 <- cbind(nmds1.2021.1.df, nmds2.2021.1) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
# ggplot(nmds_2021_dat.1, aes(NMDS1, NMDS2, color = grazing_hist, shape = nut_trt)) +
# geom_point(size = 4, stroke = 0.5) +
# labs(title = "2021 - Grouping by Grazing History & Nutrient Treatment")
ggplot(nmds_2021_dat.1, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (LOW) - Grouping by Nutrient Treatment")
ggplot(nmds_2021_dat.1, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (LOW) - Grouping by Precipitation Treatment")
Stress
0.08065034
# relativize data with Bray-Curtis dissimilarity
matrix.2021.2 <- matrix0 %>%
filter(yr == 2021,
grazing_hist == "high")
row.2021.2 <- paste(matrix.2021.2$nut_trt, matrix.2021.2$ppt_trt, matrix.2021.2$yr, matrix.2021.2$grazing_hist, sep = "_")
matrix.2021.2 <- matrix.2021.2[, c(5:79)]
rownames(matrix.2021.2) <- row.2021.2
matrix.bray.2021.2 <- decostand(matrix.2021.2, "total")
nmds.2021.2 <- metaMDS(matrix.bray.2021.2, k=2, trymax = 25)
nmds.2021.2
stressplot(nmds.2021.2)
plot(nmds.2021.2, type="t")
dummy <- as.data.frame(scores(nmds.2021.2, choices=c(1), display=c("sites")))
nmds1.2021.2.df <- as.data.frame(scores(nmds.2021.2, choices=c(1), display=c("sites"))) %>%
mutate(trt = rownames(dummy)) %>%
separate(trt, c("nut_trt", "ppt_trt", "pr", "grazing_hist"), sep="_")
nmds2.2021.2 <- as.data.frame(scores(nmds.2021.2, choices=c(2), display=c("sites")))
nmds_2021_dat.2 <- cbind(nmds1.2021.2.df, nmds2.2021.2) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_2021_dat.2, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (HIGH) - Grouping by Nutrient Treatment")
ggplot(nmds_2021_dat.2, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (HIGH) - Grouping by Precipitation Treatment")
Stress: fair
0.1230831
nmds.2021 <- metaMDS(matrix.bray.2021, k=2, trymax = 25)
nmds.2021
stressplot(nmds.2021)
plot(nmds.2021, type="t")
nmds1.2021 <- as.data.frame(scores(nmds.2021, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 3),
grazing_hist = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 4)
)
nmds2.2021 <- as.data.frame(scores(nmds.2021, choices=c(2), display=c("sites")))
nmds_2021_dat <- cbind(nmds1.2021, nmds2.2021) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_2021_dat, aes(NMDS1, NMDS2, color = grazing_hist, shape = nut_trt)) +
geom_point(size = 4, stroke = 0.5) +
labs(title = "2021 - Grouping by Grazing History & Nutrient Treatment")
ggplot(nmds_2021_dat, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 - Grouping by Nutrient Treatment")
ggplot(nmds_2021_dat, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 - Grouping by Precipitation Treatment")
# relativize data with Bray-Curtis dissimilarity
matrix.2021.block <- matrix0.block %>%
filter(yr == 2021)
row.2021.block <- paste(matrix.2021.block$nut_trt, matrix.2021.block$ppt_trt, matrix.2021.block$yr, matrix.2021.block$block, sep = "_")
matrix.2021.block <- matrix.2021.block[, c(5:79)]
rownames(matrix.2021.block) <- row.2021.block
matrix.bray.block.2021 <- decostand(matrix.2021.block, "total")
Stress: fair
0.1702944
nmds.2021.block <- metaMDS(matrix.bray.block.2021, k=2, trymax = 25)
nmds.2021.block
stressplot(nmds.2021.block)
plot(nmds.2021.block, type="t")
nmds1.2021.block <- as.data.frame(scores(nmds.2021.block, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.2021.block), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.2021.block), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.2021.block), "_"), `[`, 3),
block = sapply(strsplit(rownames(matrix.2021.block), "_"), `[`, 4)
)
nmds2.2021.block <- as.data.frame(scores(nmds.2021.block, choices=c(2), display=c("sites")))
nmds_2021_dat.block <- cbind(nmds1.2021.block, nmds2.2021.block) %>%
select(nut_trt:block, NMDS1, NMDS2)
ggplot(nmds_2021_dat.block, aes(NMDS1, NMDS2, color = block, shape = nut_trt)) +
geom_point(size = 4, stroke = 0.5) +
labs(title = "2021 (Block) - Grouping by Block & Nutrient Treatment")
ggplot(nmds_2021_dat.block, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block) - Grouping by Nutrient Treatment")
ggplot(nmds_2021_dat.block, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block) - Grouping by Precipitation Treatment")
# relativize data with Bray-Curtis dissimilarity
matrix.2021.block.low <- matrix0.block %>%
filter(yr == 2021,
block == 3 | block == 4)
row.2021.block.low <- paste(matrix.2021.block.low$nut_trt, matrix.2021.block.low$ppt_trt, matrix.2021.block.low$yr, matrix.2021.block.low$block, sep = "_")
matrix.2021.block.low <- matrix.2021.block.low[, -c(0:4)]
rownames(matrix.2021.block.low) <- row.2021.block.low
matrix.bray.block.2021.low <- decostand(matrix.2021.block.low, "total")
Stress: fair
0.1282296
nmds.2021.block.low <- metaMDS(matrix.bray.block.2021.low, k=2, trymax = 25)
nmds.2021.block.low
stressplot(nmds.2021.block.low)
plot(nmds.2021.block.low, type="t")
nmds1.2021.block.low <- as.data.frame(scores(nmds.2021.block.low, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 3),
block = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 4)
)
nmds2.2021.block.low <- as.data.frame(scores(nmds.2021.block.low, choices=c(2), display=c("sites")))
nmds_2021_dat.block.low <- cbind(nmds1.2021.block.low, nmds2.2021.block.low) %>%
select(nut_trt:block, NMDS1, NMDS2)
ggplot(nmds_2021_dat.block.low, aes(NMDS1, NMDS2, color = block, shape = nut_trt)) +
geom_point(size = 4, stroke = 0.5) +
labs(title = "2021 (Block; LOW) - Grouping by Block & Nutrient Treatment")
ggplot(nmds_2021_dat.block.low, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block; LOW) - Grouping by Nutrient Treatment")
ggplot(nmds_2021_dat.block.low, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block; LOW) - Grouping by Precipitation Treatment")
# relativize data with Bray-Curtis dissimilarity
matrix.2021.block.high <- matrix0.block %>%
filter(yr == 2021,
block == 1 | block == 2)
row.2021.block.high <- paste(matrix.2021.block.high$nut_trt, matrix.2021.block.high$ppt_trt, matrix.2021.block.high$yr, matrix.2021.block.high$block, sep = "_")
matrix.2021.block.high <- matrix.2021.block.high[, -c(0:4)]
rownames(matrix.2021.block.high) <- row.2021.block.high
matrix.bray.block.2021.high <- decostand(matrix.2021.block.high, "total")
Stress: fair
0.1412997
nmds.2021.block.high <- metaMDS(matrix.bray.block.2021.high, k=2, trymax = 25)
nmds.2021.block.high
stressplot(nmds.2021.block.high)
plot(nmds.2021.block.high, type="t")
nmds1.2021.block.high <- as.data.frame(scores(nmds.2021.block.high, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 3),
block = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 4)
)
nmds2.2021.block.high <- as.data.frame(scores(nmds.2021.block.high, choices=c(2), display=c("sites")))
nmds_2021_dat.block.high <- cbind(nmds1.2021.block.high, nmds2.2021.block.high) %>%
select(nut_trt:block, NMDS1, NMDS2)
ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2, color = block, shape = nut_trt)) +
geom_point(size = 4, stroke = 0.5) +
labs(title = "2021 (Block; HIGH) - Grouping by Block & Nutrient Treatment")
ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block; HIGH) - Grouping by Nutrient Treatment")
ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block; HIGH) - Grouping by Precipitation Treatment")
Greenhouse traits
Field traits
“H1: Pairing field results with a greenhouse 15 N experiment will allow us to relate species responses to their ability to uptake organic N”
“In addition, we will conduct a greenhouse experiment to isolate the mechanism (mineralized N vs organic N uptake rates ) by which compost may affect species composition”
Workflow plan:
overlay species trait data as vectors on the NMDS
only greenhouse traits available (only sla data from field)
use PC scores as composite trait axes OR
average the traits for each species
rda(traits.gh, scale = TRUE)
Error in colMeans(x, na.rm = TRUE) : 'x' must be numeric
Trait data is identical. Will clean data in
# dropbox.main <- read.csv("https://www.dropbox.com/scl/fi/c8vk31a807xg48jez1j03/GreenhouseTraits.csv?rlkey=ca736ea27k79mo7h8dbn67qiq&st=0h1bs53j&dl=1")
#
# dropbox.gh <- read.csv("https://www.dropbox.com/scl/fi/67aaubdydiy4ixffmwvfe/GreenhouseTraits.csv?rlkey=rm2eqeqome7hgoedhnal0gx6d&st=vkym6evx&dl=1")
#
# all.equal(dropbox.main, dropbox.gh)
# identical(dropbox.main, dropbox.gh)
# make trait species ID match matrix data (notes shown below) -------
traits.cleaned <- traits.cleaned %>%
mutate(
ID = ifelse(
ID == "Agoseris", "AGSP", ifelse(
ID == "AVEBAR", "AVBA", ifelse(
ID == "BRDR", "BRDI", ID
)
)
)
)
# summary(traits.cleaned)
traits.means <- traits.cleaned %>%
group_by(ID) %>%
summarise(
across(4:(ncol(traits.gh) - 1), list(mean = ~mean(.x, na.rm = TRUE)
)))
# make site by species matrix match traits species
matrix.names <- data.frame(colnames(matrix.bray)) #species comp (75 count)
trait.names <- data.frame(traits.means$ID) #trait species (80 count)
intersect.names <- intersect(colnames(matrix.bray), traits.means$ID) # same species (43 count)
setdiff(colnames(matrix.bray), traits.means$ID)
[1] "ASSP" "AST1" "AST2" "AST3" "CIQU" "CRTI" "CYEC" "DACA" "DIMU" "GAGR" "GAPA" "GAPH" "HYSP" "JUBU" "LYHY"
[16] "MAD1" "MAD2" "NAPU" "PAVI" "SABI1" "SABI2" "SIGA" "THGR" "TRSP" "UNBU" "UNF1" "UNF3" "UNF4" "UNF8" "UNFR"
[31] "UNGR" "XAST"
setdiff(traits.means$ID, colnames(matrix.bray))
[1] "ACMI" "ANARp" "AVFA" "BRCA" "BRHOp" "BRNIf" "BRNIp" "CLPUf" "CLPUp" "CYDA"
[11] "Crassula" "ESCA" "FEMI" "FEMY" "GITRf" "GITRp" "HYGL" "LACA" "LENIf" "LENIp"
[21] "LOPU" "LUBI" "Leontodon" "MAELf" "MAELp" "MICA" "NEMA" "PLERf" "PLERp" "RUPU"
[31] "Sanicula" "TRCI" "TRHIpi" "TRWIf" "TRWIp" "TRWIpi" "VIVA"
# check matrix codes to see if same species are named different things
unique.matrix.species <- cover.dat %>%
dplyr::select(code4, species) %>%
summarise(
code4 = unique(code4),
species = unique(species)
)
Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in dplyr 1.1.0.
Please use `reframe()` instead.
When switching from `summarise()` to `reframe()`, remember that `reframe()` always returns an ungrouped data frame and adjust accordingly.
unique.trait.species <- traits.cleaned %>%
dplyr::select(Taxon, ID) %>%
distinct(ID, .keep_all = TRUE)
change Trait name: ‘Agoseris’ –> ‘Agoseris unknown’; assuming the same as ‘AGSP’ –> ‘Agoseris sp.’ in matrix data (change to matrix name)
What’s the difference between:
‘ANAR’ and ‘ANARp’ in trait data?
‘BRHO’ and ‘BRHOp’ in trait data?
‘BRNIf’ and ‘BRNIp’ in trait data? (not in matrix data)
‘CLPUf’ and ‘CLPUp’ in trait data? (not in matrix data)
‘GITRf’ and ‘GITRp’ in trait data? (not in matrix data)
‘LENIf’ and ‘LENIp’ in trait data? (not in matrix data)
‘TRHI’ and ‘TRHIpi’ (same sp name) in trait data? (TRHI in matrix data)
‘TRWIf’ and ‘TRWIp’ and ‘TRWIpi’ in trait data? (not in matrix data)
‘Crassula’ and ‘Crassula tillaea’ (could they be the same?)
‘Sanicula’ (unknown) and ‘SABI1’ or ‘SABI2’ (could they be the same?)
‘MAELf’ (trait data) and ‘MAD1’ (Madia sp. 1 (“tall tarweed”)’ (could they be the same?)
‘MAELp’ (trait data) and ‘MAD2’ (Madia sp. 1 (“small tarweed”)’ (could they be the same?)
‘Trifolium eriocephalum’ (TRER) and ‘Triphysaria eriantha’ (TRER) (could they be the same?)
change ‘AVEBAR’ (trait data) to ‘AVBA’ (matrix data)
change BRDR –> BRDI
Stress:
0.2011127
# trait vectors onto the NMDS
envfit(nmds.comp.trait, traits.means.comp, permutations = 999)
Warning: longer object length is not a multiple of shorter object lengthError in vectorfit(X, P, permutations, strata, choices, w = w, ...) :
input data have non-matching numbers of observations
ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, fill = nut_trt, color = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Nutrient Treatment") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, color = as.factor(yr), fill = as.factor(yr))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Year")
ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, fill = as.factor(ppt_trt), color = as.factor(ppt_trt))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Precipitation Treatment") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, fill = as.factor(grazing_hist), color = as.factor(grazing_hist))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Grazing History")
ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, color = as.factor(grazing_hist), fill = as.factor(yr))) +
geom_point(aes(fill = as.factor(yr)), shape = 21, size = 4, stroke = 1) +
stat_ellipse(aes(color = as.factor(yr)),
geom = "polygon",
fill = NA,
size = 1) +
labs(title = "Grouped by Grazing History & Year") +
scale_color_manual(
values = c(
"low" = "cyan",
"high" = "brown",
"2019" = "salmon",
"2020" = "limegreen",
"2021" = "cornflowerblue")
)
ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, color = nut_trt, fill = as.factor(yr))) +
geom_point(shape = 21, size = 4, stroke = 1) +
stat_ellipse(aes(color = as.factor(nut_trt)),
geom = "polygon",
fill = NA, # Makes the ellipse transparent inside
size = 1) +
labs(title = "Grouped by Nutrient Treatment & Year") +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
# setdiff(colnames(matrix.bray), traits.gh.means$ID)
# setdiff(traits.gh.means$ID, colnames(matrix.bray))
traits.gh.means2 <- traits.gh.means %>%
filter(ID %in% intersect(colnames(matrix.bray), traits.gh.means$ID))
rownames(traits.gh.means2) <- as.character(traits.gh.means2$ID)
matrix.bray %>%
select(all_of(intersect(colnames(matrix.bray), traits.gh.means$ID))) # is this appropriate? should i relativize it all only with these species instead of just cutting out species?
nmds.comp.2 <- metaMDS(matrix.bray %>%
select(all_of(intersect(colnames(matrix.bray), traits.gh.means$ID))), k=2, trymax = 25)
nmds.comp.2 # pretty not great stress
stressplot(nmds.comp.2)
plot(nmds.comp.2, type="t")
envfit(nmds.comp.2, traits.gh.means2, permutations = 999)